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Abstract. In this paper, we will first introduce a novel multiscale represen- 
tation (MSR) for shapes. Based on the MSR, we will then design a surface 
inpainting algorithm to recover 3D geometry of blood vessels. Because of the 
nature of irregular morphology in vessels and organs, both phantom and real 
inpainting scenarios were tested using our new algorithm. Successful vessel 
recoveries are demonstrated with numerical estimation of the degree of arte- 
riosclerosis and vessel occlusion. 

1. Introduction 

1.1. Literature Reviews and Motivations. Multiscale representation (MSR) 
of functions, e.g. wavelets, has been extensively studied in the past twenty years 
[1, 2]. However, when one deals with shapes, e.g. biological shapes in R 3 , most of the 
classical theories and algorithms cannot be directly extended. In this paper, we will 
propose a new MSR for shapes based on PDEs and level set method. Although we 
shall focus on studying 3D biological shapes/surfaces, the MSR that we introduce 
here applies to general shapes/surfaces in both 2D and 3D. 

Many attempts have been made in the past on designing wavelet-typed MSR 
for 3D shapes [3, 4, 5, 6]. Among them, the method proposed by Nain et. al. 
[3, 4] is especially effective to study biological shapes. They first map the shape 
(triangulated) onto the unit sphere so that one obtains a vector-valued function 
/ : § 2 R 3 ; then apply spherical wavelet decomposition [8] to each component of 
/. However, the wavelet coefficients are not intrinsic to the shape, but dependent 
on the mapping /. Furthermore, finding a good mapping from a shape to the unit 
sphere (or to some other canonical domains) is nontrivial and in fact a popular 
ongoing research area (see e.g. [12, 13, 14, 15, 16, 17, 18, 19, 20, 21]). 

Another interesting approach was proposed by Pauly et. al. [7], where they 
introduced an MSR for point-based surfaces. Their idea was to use Moving Least 
Square method [22] to define a series of smoother and smoother point-based sur- 
faces, and then define wavelet coefficients as the displacements from two successive 
levels. Their method only requires a local parametrization of the point-based sur- 
face which is easy to calculate. However, the application of their method is rather 
limited in medical image analysis, because most of the biological shapes are not 
point-based. 
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Motivated by Pauly et. al.'s work, we will propose a new MSR for shapes in 
Section 2. The basic idea is using level set motions via solving some properly cho- 
sen Hamilton-Jacobi (HJ) like equation to obtain a sequence of shapes that become 
smoother and smoother as time evolves (analogous to coarse level approximation in 
wavelet decomposition) . Then we carefully define the so-called "details" (analogous 
to wavelet coefficients) of the MSR which carry important geometric information 
and facilitate a perfect reconstruction. While the wavelet based multi scale decom- 
position and reconstruction use filters, which are linear processes, the proposed new 
MSR for shapes uses (nonlinear) PDEs for both decomposition and reconstruction. 
However, the spirit is the same, i.e. separate features from smooth components of 
the surface. Due to the level set formulation, parametrization is no longer needed. 

1.2. Shape Modelling and Evolution PDEs. Throughout this paper, shapes 
are defined to be smooth boundaries of domains Slsl 3 and are represented by level 
set functions, typically signed distance functions. We note, however, that point- 
based and triangulated surfaces can also be handled in a similar way, as noted in 
(4) of Remark 2.2. 

A level set function <j> that represents the shape dVl is defined as follows 



We always assume that the function <f> is at least Lipschitz continuous. 

Level set motions can be achieved by solving the following H J like equation [23] , 



where we take (x,t) E V x [0, T] with V some bounded domain in R 3 and T > 0. 
Here v n (V^>) is the normal velocity, which essentially depends on while second 
order derivatives of <j) may be involved (e.g. mean curvature). If it only depends 
on first order derivatives, then it is a HJ equation. We also assume that the PDE 

(1.1) is geometric [32, 33], which guarantees contrasts invariance. Comprehensive 
theoretical analysis of PDE (1.1) and surface evolution equations can be found in 
[32, 33, 28, 29, 30, 31, 25, 26, 34]. 

The choice of velocity fields is very important and yet very non-unique. We need 
to choose one that generates a "meaningful" MSR , e.g. a sparse representation, 
for a given smooth shape. Generally speaking, we want the zero level set of u(x, t) 
becomes smoother and smoother as t increases. This is in fact a typical scale space 
behavior that has been studied for decades (see e.g. [35, 36]). It is known [35, 36] 
that under some general axiomatic hypothesis and some invariance (i.e. rotation 
and contrast invariance) assumptions on {u(x, i)}t>o, u(x,t) must be a viscosity 
solution to a PDE of the form (1.1), with the velocity field v n only depending on the 
principle curvatures of level sets of u and time t. In other words, a "meaningful" 
velocity field must be curvature dependent. 

The velocity fields that we shall focus in this paper are 

(1.2) v n = c — Xk, A > and v n = n a — k, 

where ceRis some constant, k is the mean curvature defined as k :— V • y^jj an d 
n a is the average mean curvature [38]. Note that for v n = n a — k, the PDE (1.1) 
generates an volume preserving mean curvature motion [38, 39, 40, 41]. 




x e Ct ; 

x e n c . 



(1.1) 



4> t + V n (V0) |V^| = 0, 4>{x,Q) = (j> {x) 
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2. Level Set Based MSR of Shapes: Continuous Transforms and 

Discrete Algorithms 

Let il t £ M 3 be some domain with scale t, and St ■= dVl t be the shape at scale 
t represented by some time-dependent level set function <j)(x,t), i.e. <f>(x,t) < for 
x € fit, <f>{x, t) > for x € Qf, and 

(2.1) S t = {x EM. 3 | 4>{x,t)=0} t > . 

Here So denotes the original shape with the corresponding level set function <j>o(x) = 
4>(x,0). Throughout the rest of the paper, the function <f>(x,t) is always taken to 
be the solution of (1.1). For some properly chosen v n in (1.1), e.g. with v n = —k 
or n a — k, we can obtain a continuous series of shapes {St}te[o,T], which tends to 
become smoother when t increases. Based on this, we define our continuous level 
set based MSR of So as follows. 

Definition 2.1. Let<j>(x,t) be the solution of the PDE (1.1) and (x,t) G Vx [0,T]. 
We now understand xi(t) as a path on the propagating l-th level set of <f), i.e. 
(f>(xi(t),t) = I. For simplicity, we shall omit the subscript "I" unless a particular 
level set is considered. 

(1) We now define the multiscale transformation (MST) of <po(x) as 

(2.2) W(x, t) := W{4>o) := "«» = -A*)- 

Vector —x'(t) is the displacement vector and w(x,t) := —v n (x,t) is the 
detail of the MST. 

(2) We shall call W(x, t) the displacement vector field at scale t, and denote 
W\{x,t) (u>i(x,t)) as the restriction ofW{x,t) (w(x,t)) on S t . 

(3) The MSR for the original shape S in terms of (j) (x) is denoted as 

MSR(4>o) = {{W{x,t)} tmT) ,(j>{x,T)). 

(4) We define the inverse multiscale transformation (IMST) via solving 
the following PDE 

(2.3) Vv + VF(x, T - r) • VV> = 0, i/)(x,0) = 4>{x,T). 
for given T > and < r < T. 

Remark 2.2. 

(1) The technique of generating a sequence of the spaces {St} via solving PDEs 
is known as scale space decomposition (see e.g. [35, 36]j. However, a 
classical scale space analysis does not study the details as defined in (2), 
and does not have a reconstruction as in (4). 

(2) The last identity in (2.2) can be easily shown by using PDE (1.1) and the 
assumption that x'(t) is aligned with normal directions of level sets of <p. 

(3) The detail w\{x, t) is a function on St that characterizes intrinsic geometric 
information of the shape at scale t. Here by intrinsic we mean that w\ [x, t), 
as well as {St}t>a, does not depend on the initial embedding <fio for a large 
class of functions [33], but only depends on So- Therefore, we now have an 
intrinsic MSR for So •' 

(2.4) MSR(S ) - {{W>|(x,i)} te(0 , T) ,S T }. 
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Furthermore, the above MSR is invariant under translation and rotation of 
So. 

(4) The MSR defined above can be easily adapted to a point-based or triangu- 
lated surface. One simply need to first associate the surface with a level set 
function and then perform the MST. For point-based surfaces, the IMST 
from its MSR (2.4) can be point-wise defined as So = St + J W\(x,i)dt 
or equivalently xo(0) = xq{T) + —x' Q (t)dt, which is obviously true. 

Now the question is that if we have perfect reconstructions via (2.3). The answer 
is given in the following proposition, which directly follows from theories of ODEs. 

Proposition 2.3. Assume that W(x, t) stays Lipschitz continuous for {x, t) € T> x 
[0, T]. Then the equation (2.3) inverts the MST defined by (2.2) in the sense that 
ip(x,r) :— <j)(x,T — t) is the unique solution of (2.3). 

Remark 2.4. 

(1) The assumption in Proposition 2.3 is not always valid (e.g. v n = c < and 
</>o(x) representing a cube). However, if we choose, for example, v n = —n, 
v n = K a — k or v n = c — \k (for A > 0), and choose some appropriate 
ending time T > (e.g. before any topological changes occur), the above 
assumption will be valid and we will have a perfect reconstruction using 
(2.3) [33, 37]. 

(2) Generally speaking, the vector field W(x,t) does not stay Lipschitz globally 
in time. This happens precisely when the corresponding surface evolution 
starts to develop singularities, and it is difficult to find a surface evolution 
that guarantees to have global smooth solutions for a general initial surface 
So- For some special class of initial surfaces, however, it is relatively easy to 
find such motion. Taking v n = n a — k for example, it is shown in [38] that if 
the initial surface So is close enough to a certain sphere (but not necessarily 
convex), then St stay smooth and eventually converges exponentially fast to 
a sphere with the same volume as So- 

Notice from Definition 2.1 and Proposition 2.3 that to perfectly reconstruct 4>o(x) 
from cj)(x, T), we need to store the entire vector field W(x, t) for every ieP and all 
scale t. However, in practice, we only want a perfect reconstruction of So, and thus 
we do not need that much information. Therefore, only the displacement vectors 
within a narrow band of the zero level set of <fi(x, t) need to be stored. 

We can be even more "greedy" here by only storing W\(x,t). When performing 

inverse transform, we will need to extend W\ (x, t) to at least a narrow band of the 
zero level set of (p(x,t). Note that no extension can guarantee an exact recovery 
of the vector field W(x,t), and hence the reconstruction of So will not be exact. 
However, if the extension is conducted accurately and the mesh grid is dense enough, 
i.e. the resolution of the shape is high enough, the reconstruction should be more 
and more accurate. The extension we shall adopt here is such that the extended 
vectors are constant in the normal directions of each level set of (j)(x,t) [42]. For 
simplicity, we will use a local search method to extend W\ (x, t) to a narrow band 
of the zero level set of (f)(x, t). 

Our proposed discrete version of MSR is given in Algorithm 1. 
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Algorithm 1 Level Set Based MST and IMST 

Start from the given level set function (j>o(x) representing shape So- Choose time 
steps = to < ti < . . . < t/v = T, where maxj(tj + i — ti) is small. 
Initialize: Sample a point set Xq from So (cither uniformly or non- uniformly) . 

MST: 

while i < N do 

1. Starting from 4>(x, solve PDE (1.1) for t e [tj_i, tj] and obtain </>(#, ti). 

2. Orthogonally project onto the zero level set of (j){x,ti) and obtain Aj. 
2. Compute the discrete displacement vector by W\i = X^— Aj_i, and i <— 

end while 

We then obtain the discrete MSR of S : MSR(Sb) := 
{W ll ,W l2 ,...,W ]N ^(x,T)}. 

IMST: 

1. Extend the vector fields {W^}^ such that the values are constant along 
normal directions of the level sets of <p(x, ti). 

2. Solve (2.3) using W\i within interval [tj, ti— i] iteratively for each i. 



3. Numerical Experiments on the MSR 

One of the key steps of implementing Algorithm 1 is to solve the evolution 
PDE (1.1) efficiently. There are many ways of solving equation (1.1). The most 
straightforward way is to use monotone finite difference schemes [23, 24]. However, 
it is not very efficient computationally. To overcome this, Merriman, Bcnce and 
Osher introduced a diffusion-based level set motion in [43, 44], and it was further 
studied in [45, 46, 47, 48], where in [45] the correctness of the method is rigorously 
proven. In [40], Ruuth and Wetton introduced a fast algorithm to calculate volume 
preserving motion by mean curvatures. All these methods speeded up curvature 
driven motions drastically. 

In this section, we will recall the fast algorithms of level set motion for the case 
v n = c and v n = n a — n given by [40, 43, 44, 47]. These algorithms will be used in 
the later sections to generate fast multiscle decompositions of shapes. 

Now we first recall the fast method of solving (1.1) with v n — c (see [43, 44, 47]) 
in Algorithm 2. 



Algorithm 2 Level Set Motion with Constant Normal Velocity 

Start from a given shape represented by <j>. 
while t < T do 

1. Define the corresponding characteristic function by \ — 1{^ <0 }. Set Vq 
equal to the volume of {<j> < 0}. 

2. Starting from \, evolve \ f° r a time At by \t — V 2 %- 

4. Sharpen: 

f 1 ifx>0 
[0 otherwise 

5. Let t <— t+ At. Compute <p(x,t) from \ y i a fast sweeping method [49]. 
end while 
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We now recall the fast implementation of (1.1) with v n = n a — k proposed by 
Ruuth and Wetton [40] in Algorithm 3. Their algorithm is based on the diffusion- 
based mean curvature motion proposed by [43, 44] . Note that if we remove step 3 
in Algorithm 3 and choose A = 0.5 in step 4, it is exact the fast mean curvature 
motion proposed in [43, 44]. 



Algorithm 3 Volume Preserving Mean Curvature Motion: v n = n a — n. 

Start from a given shape represented by 4>. 
while t < T do 

1. Define the corresponding characteristic function by \ — 1{0 <O }- Set Vo 
equal to the volume of {0 < 0}. 

2. Starting from y, evolve x f° r a time At by Xt = V 2 %- 

3. Determine the threshold value that preserves the volume of the set: i.e. 
find a < A < 1 s.t. 



| otherwise 

5. Let t <— t + At. Compute <p(x,t) from x via fast sweeping method [49]. 
end while 



Some numerical results of the MST and IMST in Algorithm 1 are presented in 
Figure 1 using the tested shape (right hemisphere of a cortex). The velocity field 
is chosen to be v n = n a — n and 5 levels of decomposition are conducted (first and 
second row of Figure 1). Details Wu are drawn on the surface Si (second row of 
Figure 1), where the value is positive, when Wu is pointing outwards and negative 
when it is pointing inwards. The IMST is also presented in Figure 1 where Si 
denotes the reconstruction of level i from level i + As we can see, although the 
reconstructions are not exact for each level, they are quite accurate in the sense 
that most of the features are well reconstructed. We also illustrate sparseness of 
the coefficients {VF|i}f =1 in Figure 2, where one can see that the energy of W\i are 
concentrated around 0, especially for the later levels. 



Evaluating missing parts in medical images provides important information as 
signs of diseases. One of the most common situation is the phenomenon of vessel 
narrowing or occlusion in angiographic images. Estimating these abnormalities can 
help document disease progression. 

The recovery of image surfaces such as blood vessels can be regarded as a surface 
inpainting problem [64, 65]. Inpainting problems, for both images and surfaces, have 
been extensively studied in the literature [50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 
61, 62, 63, 64, 65, 66, 67, 68, 69, 70]. They occur when part of the data in an image 
or regions of a surface is missing or corrupted. The major task of inpainting is to 
fill in the missing information based on the geometry of the image/surface. In this 
section, we will propose a new surface inpainting algorithm based on the MSR in 
Section 2 for blood vessel reconstruction that arises in medical image analysis. 



\{x: X < A}| - Vb <e. 



4. Sharpen: 




4. Application in Blood Vessel Recovery 
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FIGURE 1. First row (left to right): MST S , S 1: . . . , S 5 . Second 
row shows the details of MSR on Si, . . . , S5. Third row shows 
IMST Si, i = 0,1,..., 4, where the Hausdorff distance between 
Si and S*, are: l.Uh, 0.7ih, 0.74ft,, 0.69ft, and 0.63ft respectively 
(with ft the mesh size). 
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FIGURE 2. Histograms of Wu for i = 1, . . . , 5 (left to right). 



Our surface inpainting algorithm (Algorithm 4 below) inherits the following 
framelet-based image inpainting structure proposed by Cai et. al. [61]: 

(1) Take framelet transform of the given image; 

(2) Truncate the framelet coefficients via soft-thresholding and reconstruct; 

(3) Apply the exact data outside the inpainting domains, and repeat. 

Since we already have an MSR for surfaces, the first step above can be replace by 
our MST. For the second step, we shall solve the following PDE for IMST instead 
of the PDE (2.3) that was originally proposed in Definition 2.1: 

(4.1) ip T + W(x, T - t) ■ VV> = eVV, ^(a;, 0) = <f>(x, T). 

The above PDE mimics thresholding in the sense that it penalizes the reconstruction 
from W by introducing a vanishing viscosity eV 2- 0, which forces some information 
outside the inpainting region flows into the inpainting regions. Also, when e — > 0, 
the solution of (4.1) converges to the viscosity solution of (2.3) [26, 34]. 
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Since we generally expect volumes of surfaces to increase during inpainting, we 
choose the following PDE for the MST, 

(4.2) 4> t + {c + K a - K)\V<j)\ =0, <p(x,0) = <f>o(x), c >0- 

Note that the PDE (4.2) generates a mean curvature motion with increasing vol- 
umes of the domains enclosed by level sets of <j>(x, t). The constant c can be regarded 
as a parameter that needs to be adjusted according to different surface inpainting 
scenarios. In our experiments, we solve PDE (4.2) via a combination of Algorithm 
2 and Algorithm 3 recalled in Section 3 



Algorithm 4 Surface Inpainting via MSR 

Start from <po, with inpainting region D. Choose some e > 0. 
while "Not converge" do 

1. Perform discrete MST by solving (4.2) and acquire W\i by Algorithm 1. 

2. Perform IMST by solving (4.1) and obtain tp s . 

3. Copy the known information to ip e : i\} e \D<^ iPo\d<=- 

4. Decrease amount of smoothing: e \. 
end while 



We test Algorithm 4 on both phantom (first two vessels in Figure 3) and real 
(last two vessels in Figure 3) surface inpainting scenarios. First row of Figure 3 
shows four blood vessels with inpainting regions specified by red circles. For the 
two phantom inpainting scenarios, the inpainting regions are created manually, and 
the surface within those regions were chopped off. For the two real inpainting 
scenarios, we do not know the exact inpainting regions. Therefore in practice, we 
adopt a user interactive strategy to determine the inpainting regions. After several 
points have been selected on the surface, the inpainting regions are then generated 
automatically. Inpainting results are given in second and third row of Figure 3. 
We want to point out that during the inpainting process, topological change may 
occur for some cases (e.g. second vessel in Figure 3). Although it violates the 
assumption in Proposition 2.3, topological change is still allowed for inpainting 
problems. The reason is that perfect reconstruction is only required at the very 
last stage of inpainting (i.e. when e w 0) in order to ensure convergence, while 
topological changes will happen during the middle of the process if the parameters 
(e.g. c in (4.2)) are properly chosen. 

5. Conclusion 

In this paper, we introduced a novel multiscale representation (MSR) for shapes 
which is intrinsic to the shape itself, does not need any parametrization, and the 
details of the MSR reveals important geometric information. Based on the MSR, 
we then proposed a surface inpainting algorithm and applied it to recover corrupted 
blood vessels. This technique is especially useful to study arteriosclerosis and vessel 
occlusions. Numerical results showed that the inpainting regions were nicely filled 
in according to the neighboring geometry of the vessels and allowed us to accurately 
estimate the volume loss of vessels. 

There are still many interesting aspects of both the MSR itself and its appli- 
cations worth discovering. For example, a rigorous analysis of how Algorithm 1 
approximates the continuous version in Definition 2.1 needs to be done. Also, 
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FIGURE 3. Blood vessel inpainting. Row 1: vessels before inpaint- 
ing; row 2: vessels after inpainting; row 3: inpainted regions shown 
in red. The percentages of the volume of inpainted region over that 
of the entire shape are: 5.3%, 19.2%, 6.7% and 5.7%. 

we can apply our MSR framework to other type of shape processing and analysis 
problems, e.g. shape registration and classification. In fact we believe that many 
techniques based on classical MSR (or multiresolution analysis) for functions (e.g. 
wavelets) can now be extended to shapes. 
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